double gaussrand(void)
{
	static double U, V;
	static int phase = 0;
	double Z;
	
	if(phase == 0)
	{
		U = rand() / (RAND_MAX + 1.0);
		V = rand() / (RAND_MAX + 1.0);
		Z = sqrt(-2.0 * log(U))* sin(2.0 * M_PI * V);
	}
	else
	{
	Z = sqrt(-2.0 * log(U)) * cos(2.0 * M_PI * V);
	}
 
	phase = 1 - phase;
    return Z;
}
